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FINITE  DIFFERENCE  METHODS  FOR  TIME  DEPENDENT, 
LINEAR  DIFFERENTIAL  ALGEBRAIC  EQUATIONS' 


BY 

Patrick  J.  Rabier  and  Werner  C.  Rheinboldt^ 


Abstract.  Recently  the  authors  developed  a  global  reduction  procedure  for  linear,  time- 
dependent  DAE  that  transforms  their  solutions  into  solutions  of  smaller  systems  of  ODE’s. 
Here  it  is  shown  that  this  reduction  allows  for  the  construction  of  simple,  convergent  finite 
difference  schemes  for  such  equations. 


1.  Introduction:  Time-dependent,  linear  differential  algebraic  equations  (DAEs), 

(1.1)  A{t)x  +  B{t)x  -  b{t),  A(t),D{t)  €  r(R"),  bit)  e  R",  t  €  R 

arise  in  many  circuit  and  control  problems  (see  e.g,  [C83],  [C85],  or  [CBP89]  for  some 
references).  In  general,  standard  ODE)- type,  numerical  methods  of  for  (1.1)  are  known 
to  fail  or  perform  poorly  when  the  index  exceeds  one,  and  certain  Taylor-type  methods 
developed  in  [C85],  that  apply  to  a  larger  class  of  problems,  require  the  solution  of  large- 
augmented  systems  of  equations. 

Recently  we  developed  a  new,  global  reduction  theory  ([RR93])  which  leads  to  existence 
and  uniqueness  results  for  classical  as  well  as  generalized  solutions  of  (1.1)  under  rather 
general  conditions.  The  aim  here  is  to  show  that  this  reduction  process  allows  for  the  con¬ 
struction  of  simple,  convergent  finite  difference  approximations  for  the  numerical  solution 
of  (1.1)  which  appears  to  provide  a  new  tool  for  the  solution  of  general  systems  of  the  form 

(1.1) . 

2.  Summary  of  the  Reduction  Process:  In  its  basic  form  the  reduction  process  of 
[RR93]  assumes  that  the  coefficient  functions  A  and  B  are  analytic.  Although  general¬ 
izations  to  the  non-analytic  case  are  also  discussed  in  [RR93],  we  shall  retain  here,  for 
simplicity,  the  analyticity  assiunption.  Throughout  this  note  analytic  mappings  will  be 
referred  to  as  “mappings  of  class  C'^" . 

A  main  tool  for  the  proof  of  the  globality  of  the  reduction  procedure  is  the  concept  of 
“transformation  functions”  introduced  by  T.  Kato  [K50].  With  it  and  with  another  result 
of  Kato  [K82]  the  following  basic  re.sult  was  proved: 

’The  work  was  supported  in  part  by  ONR-grant  N-0001.t-90-J-1025,  and  NSF-grant  CCR-9203488 
^Department  of  Mathematics  and  Statistics,  University  of  Pittsburgh,  Pittsburgh,  PA  15260 
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Theorem  1.  Let  J  C  R  be  an  open  interval,  m  >  1  any  integer,  and  M  €  £(R"’ )) 

with  r  =  maxtgj  rank  M{t).  Then,  there  is  a  subset  S  d  J  of  isolated  points  such  that 
rank  M{t)  =  r  if  and  only  if  t  E  J  \  S.  Furthermore,  the  orthogonal  projections  onto 
rge  M{t)  and  ker  M(t),  t  ^  jJ\S  are  analytic  on  t  ^  \  S  and  can  be  extended  as  analytic 

functions  over  the  entire  interval  J . 

With  the  notation  of  this  theorem  the  orthogonal  projection  P{t)  onto  rge  Mit).  for 
t  €  J^\S,  is  at  any  point  of  5  cin  orthogonal  projection  onto  a  subspace  containing  rge  M{t) 
and  this  subspace,  called  the  extended  range  of  M{t)  and  denoted  by  ext  rge  M{t).  is 
independent  of  the  specific  choice  of  P. 

For  ease  of  notation  we  call  a  pair  {A,D)  of  coefficent  functions  admissible  if  € 

C{W‘))  on  a  fixed  open  interval  C  R.  An  admissible  pair  is  said  to  be  regular 
if  rank  [,4(f).4(t)^  4-  B{t)D{t)'^]  =  n.^t  ^  J .  Then  Theorem  1  ensures  that  rank  A(t)  — 
r,  Vt  €  \  S  where  S  C  J  consists  of  isolated  points,  and  that  there  is  a  there  is  a 

family  of  projections  P  G  £(R”))  onto  ext  rge  A(t),  Vt  6  J-  Now  Q  =  I  —  P 

can  be  shown  to  satisfy  dim  ker  Q(t)5(t)  =  r,  'it  £  J.  Moreover,  there  exist  mappings 
C  G  C‘^(J;£(R'';R"))  and  D  G  C-'l  J; £(R",R''))  such  that 

(2.1)  C{i)  ^  GL[W, ker  Q(t)B{t)),  it  e  J ■ 

(2.2)  Z)(0  €  G£(ext  rge  A(0,R^),  it  e  J , 

respectively.  For  any  such  choice  of  mapppings  Q,  C,  D  the  pair  of  admissible  (Ai,i?]) 
defined  by 

(2.3)  Ai ,  Bi  G  C‘^(  J;  £(R")),  Aj  =  DAC,  Bj  =  D[BC  4-  AC), 

is  called  a  reduction  of  the  regular  pair  (A,B)  (in  J7).  Clearly,  such  a  reduction  is  not 
unique,  since  it  depends  upon  the  choice  of  C  and  D.  However,  any  two  reductions  of 
(A,  B)  turn  out  to  be  equivalent. 

More  specifically,  between  admissible  pairs  (A,  B)  and  (A,  B),  define  a  relation  (A.  B)  ~ 
(A.B)  in  4.7  by  the  conditions  A  =  MAN  and  B  =  M{BN  4-  AA”^)  for  some  M,N  G 
C£(R")).  This  turns  out  to  be  an  equivalence  relation  on  the  set  of  all  admissible 
pairs  (in  J).  Furthermore,  if  (A,B)  ~  (A, B)  then,  on  J .  we  have  (i)  rank  A(t)  = 
rank  A(t),  (ii)  {A,B)  is  regular  if  and  only  if  {A,B)  is  regular,  and  (iii)  If  (A,B)  and 
{A,B)  are  regular  any  reduction  (Ai,Bi)  of  {A,B)  is  equivalent  to  any  reduction  (Ai.Bj ) 
of  (A,B). 

This  permits  the  definition  of  our  reduction  procedure  which,  up  to  equivalence,  is 
independent  of  any  particular  choice  that  must  be  made  at  each  step.  Suppose  that  the 
admissible  pair  (A,  B)  is  regular.  Then,  any  two  reductions  (Aj ,  Bi )  and  ( A] ,  Bj )  of  ( A.  B) 
are  either  both  regular  or  not.  If  they  are,  then  once  again  any  reductions  (A2,B2)  of  and 
(A2,B2)  of  (Aj.Bj)  and  (Ai,B]),  respectively,  will  be  simultaneously  regular,  and  so  on. 
In  line  with  this,  an  admissible  pair  is  called  completely  regular  (in  N)  if  this  reduction 
procedure  can  be  continued  indefinitely.  In  particular,  the  complete  regularity  of  a  pair 
(A.B)  in  J  implies  its  regularity  in  J. 
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For  a  completely  regular  pair  {A,B)  consider  a  sequence  (Aj,Bj),  j  =  0.1 . 4o  =  .4. 

Bo  =  B,  such  that  {Aj.Bj)  is  some  reduction  of  (.4j_i,  ),  j  >  0.  Then.  Vj  == 

maxtgj[rank  is  independent  of  the  specific  choice  of  {Aj.Bj),  j  >  0  and  with 

r_i  =  n  we  have  Aj{t)  €  for  t  ^  J  and  j  >  0.  Moreover,  the  >  0  decrease 

monotonically  with  j  whence  there  exists  a  smallest  integer  0  <  u  <n  such  that  =  ri,_i 
and  Au{t)  G  GLiW''’-'^ ),  Vf  €  i7  \  5^  where  Si,  C  J  consists  only  of  isolated  points.  This 
integer  v  is  the  index  of  the  pair  [A,B). 

Theorem  2.  Let  (1.1)  be  a  DAE  with  an  admissible  pair  (A,B)  of  coefficient  functions 
and  any  b  G  l<k<ooork=uJ.  Suppose  that  (Ai.Bi)  is  any  reduction 

of  (A,B)  defined,  say,  by  the  projection  Q  and  the  mappings  C,  D  satisfying  (2.1 ),  (2.2). 
Then  there  exists  uq  G  C*(J7’; R”)  such  that  B{t)uo{t)  —  b{t)  G  ext  rge  A{t),  Vt  G  J, 
and  indeed  we  may  use  uq  =  B^(AA^  +  BB^)~^b.  With  this,  a  differentiable  mapping 
X  :  JT"  R"  solves  (1.1)  if  and  only  if  x  —  Cxi  +  uq  where  Xj  :  JT"  — +  R''  is  a  differentiable 
solution  of  the  reduced  DAE 

(2.4)  Ai(t)xi  +  Bi(t)x  =  6i(t),  tej,  bi  =  D{b  -  Buo  -  Aiio). 

By  recursive  application  of  this  result  it  follows  that  when  {A,  B)  is  completely  reducible 
with  index  u  then  after  n  steps  we  arrive  at  a  linear  ODE 

(2.5)  A„{t)xi,  +  By{t)x„  =  b„(i), 

where  Ai,(t)  G  GL(R’'‘’"')  except  perhaps  at  isolated  points.  Thus  (2.4)  is  equivalent  with 
an  explicit  ODE  with  singularities.  A  simple  example  for  this  is  the  system 

/  2t  2  0  \ 

I  0  0  2jx+x=0 

\  -2t^  -2t'^  -2t ) 

given  in  [CSS]  which  has  index  2.  Here  the  final  reduced  ODE  is  — 2t(l  +  3t^  +  #'^)x2  + 
(—3  —  8t^  4-  3t®)x2  =  0  and  hence  is  singular  at  t  =  0.  In  [CSS]  it  was  noted  that  the  DAE 
has  index  2  for  t  >  0  and  the  singularity  at  t  =  0  was  interpreted  as  index  3  at  that  point. 

Of  course,  the  case  of  (1.1)  is  of  special  importance  when 

(2.6)  Ai,{t)  e  ytej. 

The  condition  (2.6)  is  independent  of  the  reduction,  and  if  it  holds  then  (2.5)  is  an  explicit, 
linear  ODE  without  singularities  for  which  the  stzindard  existence  and  uniqueness  results 
are  applicable.  Thus  if  the  initial  condition  x(fo)  =  to  ^  J  satisfies  a  certain  consistency 
condition  (see  [RR93])  the  resulting  initial  value  problem  for  (1.1)  has  a  globally  defined 
unique  solution. in  JT. 

3.  Finite  Difference  Approximations:  The  global  reduction  of  the  DAE  (1.1)  allows 
for  the  construction  of  convergent  finite  difference  approximations  of  (1.1).  These  approx¬ 
imations  are  applied  to  the  reduced  DAE  at  stage  j  —  1  of  the  process  which  after  one 
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further  step  becomes  the  ODE  (2.5)  for  which  we  assume  now  specifically  that  (2.6)  holds. 
In  other  words,  the  finite  difference  scheme  is  applied  to  a  (reduced)  DAE  of  index  one. 
Thus,  for  ease  of  notation,  we  assume  here  simply  that  the  original  DAE  (1.1)  has  index 
one.  As  before  (Ai,.Bi)  denotes  a  reduction  of  the  admissible  pair  iA,B)  defined  by  the 
mappings  Q,  C,  D.  We  also  note  that  the  index  one  assumption  is  then  equivalent,  for 
instance,  with  the  invertibility  of  A{t)  +  Vt  €  J ,  (see  [RR93]). 

For  a  given,  sufficiently  small  step  h  >  0  suppose  that  t,  =  to  +  ih  6  for  ?  =  0, 1, . . .  .m 
and  consider  first  then  explicit  Euler  approximation 

(3.1)  A{ti)^{x,+i  -  a:,)  +  B{t,)x,  =  fe(t,). 

Then  the  following  solvabihty  result  holds: 

Theorem  3:  Any  solution 

(3.2)  a:o,a:i, . . .  ,Xm  €  R" 
of  (3.1)  satisfies  for  i  =  0, 1, . . .  ,  m  —  1  the  equations 

(3.3)  Q(ti)B{t,)xi  =  Q{t,)b(ti), 

[A(fi)  +  Q(tj+i)B(fi+i)]ij^.j 

(3.4)  =  [A(ti)  -  hB{ti)]xi  +  hb{ti)  +  Qiti+i  )b{ti+i ). 

Conversely,  for  sufficiently  small  h  and  any  given  xq  €  R”  satisfying  (3.3)  (at  tg),  the 
solution  (3.2)  of  (3.4)  is  unique  and  is  also  a  solution  of  (3.1 ). 

Proof:  For  any  solution  (3.2)  of  (3.1)  it  follows  after  multiplication  by  Q{ti)  that  (3.3) 
holds  for  0  <  z  <  m  —  1.  Hence,  by  adding  for  z  =  0, . . .  ,  ttz  —  1  the  (z  +  l)-st  equation  (3.3) 
to  (3.1)  we  obtain  (3.4).  The  smoothness  of  A,  B,  and  Q  ensures  that,  for  sufficiently  small 
h,  the  matrix  in  the  square  bracket  on  the  left  of  (3.4)  is  nonsingular  for  z  =  0, . . .  ,  zr?  —  1 
since,  as  noted  above,  A(t,)  +  Q{t,)B{ti)  is  invertible  by  the  index-one  assumption.  Hence 
for  given  xq  the  solution  (3.2)  of  (3.4)  is  uniquely  determined.  If  (3.3)  holds  for  xg  (at 
to)  then,  by  induction  on  z,  it  follows  that  this  solution  (3.2)  of  (3.4)  satisfies  (3.3)  for 
z  =  0, . .  .  ,  zrz  —  1.  In  fact,  if  (3.3)  is  valid  for  some  z,  0  <  z  <  m,  then  we  obtain  after 
multiplication  of  (3.4)  by  Q(^)  that  Q(t,)Q(ti+i)5(ti+i)xi+i  =  Q{ti)Q{t,+i)b(t,^i).  Since 
Q(tt)  is  injective  on  rge  (5(ti+])  =  rge  A(t,+i)'*‘  (for  t,+i  —  t,  small  enough)  this  implies 
that  (3.3)  holds  for  z  +  1  in  place  of  z.  Now  by  adding  the  z  -f  1-st  equation  (3.3)  to  (3.4) 
we  see  that  the  solution  of  (3.4)  also  solves  (3.1).  □ 

For  any  solution  of  (3.4)  it  follows  from  (3.3)  that  x,  —  ug(t,)  €  ker  Q{t,)B(t,)  for  all  z 
whence 

(3.5)  Xi  -  Ug{ti)  =  Cit,)zi,  Vz, 
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for  some  sequence  z,  €  2  =  0, . . .  ,m  +  1.  Then,  using  that  )uo(^i+i )  —  ^(^+i ) 

belongs  to  [rge  )]■'■  we  obtain  from  (3.4),  after  a  simple  calculation,  that 

A(t,)C(<.+i)^(-,+i-M)  +  [B(t.)C(f.)  +  A(f.)i(C;(t.+i)-C(t.))] 

(3.6)  =  c,  =  b(ti)  -  A(t,)^{uoit,+i)  -  uo{t,))  -  B{t,)uo{t,). 

After  multiplying  this  by  D{t,)  we  see  that  (3.6)  is  a  finite  difference  approximation  of  the 
initial  value  problem 

(3.7)  Ai{t)z  +  Bi(t)z  =  c(t),  z{to)  =  zo, 
where  zq  is  characterized  by  the  condition  C(tQ)zQ  =  xo  —  uo(fo),  and 

(3.8a)  Ai{t)  ^  D{t)A{t)Cit  +  h), 

(3.86)  BAt)  =  Dit)[B{t)C(t)  +  Ait)hc{t  +  h)  -  C(0)], 

n 

(3.8c)  b{t)  =  b{t)  -  A(<)y(uo(<  +  h)-  Wo(0)  “  B{t)uo{t). 

n 

Since  (2.6)  was  assumed  to  hold  for  the  ODE  (2.5)  obtained  at  the  last  reduction  step, 
(2.5)  is  equivalent  with  the  linear,  explicit  ODE 

(3.9)  z  =  K{t)z  +  k{t),  A'(t)  =  i,(0"'Bi(t),  k{t)  =  AAtr^b{t). 

Hence,  by  a  standard  estimate  for  Euler’s  method  (see  e.g.  [HNW87],  Theorem  7.5)  we 
obtain 

(3.10)  IUit.)-^.||  <Ki|^|- 

On  the  other  hand,  a  comparison  of  (3.7)  and  the  reduced  equation  (1.3)  shows  that  for 
6^0  the  difference  of  the  coefficient  functions  (3.8)  and  the  corresponding  coefficient 
functions  of  (2.5)  is  of  order  \h\  uniformly  on  [toi^m]-  Thus  a  standard  application  of 
Gronwall’s  inequality  provides  that 

max  llxi(#)  -  2(f)ll  <  «2|61, 

20<2<<m 

whence  altogether  we  obtain  from  (1.4;  and  (3.5)  that 

(3.11)  ||xi  -  x(t,)||  =  ||C(<,)(2ri(ti)  -  2, 011  <  k|/i|,  Vi, 

with  K  =  maxt(,<(<t„  ||C(t)||(Ki  +K2)-  In  other  words,  when  xq  satisfies  (3.3)  at  to  then  the 
unique  solution  of  the  difference  equation  (3.4)  provides  an  approximation  of  the  solution 
of  (1.1)  with  global  error  0{h). 
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The  result  is  easily  extended  to  variable  steps  in  t  by  using  a  continuous  grid  function. 
It  is  also  conceptually  strriightforward  to  derive  higher  order  discretizations  or  implicit 
schemes.  We  illustrate  this  briefly  for  the  implicit  Euler  scheme 

(3.12)  —  I,  )  +  )x,4.i  =  \  ), 

Now  any  solution  (3.2)  satisfies  for  i  =  0,l,...  ,m  —  1  the  equations 

(3.13)  Q(t,+i  )5(ti+i  )xi+i  = 

[>l(^i+i )  +  Q(^«+i  )-B(t,+i )  +  )]x,4.i 

(3.14)  =  ^(t,+i  )x,  4-  hb(ti+i )  -f  Q(t,^.i  )b(ti^i ), 

For  small  h  the  matrix  on  the  left  of  (3.14)  is  nonsingular  and  hence,  for  given  xq  the 
solution  (3.2)  of  (3.14)  is  unique.  This  solution  satisfies  (3.13)  for  i  =  0, . . .  ,  m  —  1  as  is 
readily  seen  by  multiplying  (3.14)  with  Q(<i+i)  and  dividing  by  /i  +  1.  Thus  subtraction 
of  (3.13)  from  (3.14)  shows  that  the  solution  of  (3.14)  also  solves  (3.12).  Moreover,  the 
unique  solution  of  (3.14)  represents  again  an  approximation  of  the  solution  of  (1.1)  with 
global  error  0{h).  The  proof  is  entirely  analogous  to  that  given  for  (3.4)  and  will  not  be 
repeated  here. 

If,  for  >  0  the  matrix  A(t,+i)  +  hB{ti^i)  in  (3.12)  is  nonsingular  then  the  solution  of 
(3.12)  can  be  computed  directly.  This  observation  is  well  known  to  apply  also  to  higher 
order  BDF  formulas  and  is  the  basis  of  the  widely  used  DAE  solver  DASSL  (see  [BCP89]). 
However,  the  nonsingularity  assumption  requires  the  matrix  pencil  A{t),  B{t)  to  be  regular 
for  all  t  which  is  not  necessarily  true  for  all  index  one  problems.  Moreover,  since,  in 
any  case,  the  matrix  becomes  singular  for  ft  =  0,  increasing  difficulties  are  expected  for 
decreasing  ft.  This  problem  is  not  shared  by  either  (3.4)  or  (3.14). 

Difference  schemes  as  (3.4)  and  (3.14)  open  up  surprisingly  effective  numerical  methods 
for  the  solution  of  the  systems  (1.1).  In  particular,  (3.4)  can  be  used  as  the  base  method 
in  an  explicit  extrapolation  integrator.  An  implementation  of  the  resulting  algorithm  for 
general  index-one  problems  (1.1)  has  been  called  LTVXE.  It  has  the  general  form  of  the 
extrapolation  code  EULEX  (see  [DNP88])  and,  as  EULEX,  it  is  based  on  the  order  and 
step  control  mechanism  of  [D83]. 

Of  course,  when  (1.1)  has  index  exceeding  one,  the  application  of  this  integrator  pre¬ 
supposes  the  availability  of  the  DAE  arising  at  stage  v  —  1  of  the  reduction  process.  Under 
the  assumption  that  subroutines  for  all  needed  derivatives  of  the  coefficients  A,  B,  and  b 
are  given,  a  computational  implementation  of  the  reduction  process  is  feasible.  This  will 
be  discussed  elsewhere. 

Here  we  show  only  one  simple  example  of  the  method  when  applied  to  the  index  one 
problem 
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This  is  an  example,  called  not  regular  in  [CBP89],  for  which  the  matrix  pencil  is  singular 
for  all  t  €  R.  Accordingly,  as  noted  above,  DASSL  cannot  be  expected  to  solve  (3.15),  and. 
in  fact,  DASSL  failed  consistently  for  any  tolerance  with  either  one  of  the  error  messages 
"'the  iteration  matrix  is  singular"  or  ‘Hhe  corrector  failed  to  converge  repeatedly  or  with 
abs(h)  =  hmin". 

Our  reduction  shows  readily  that  (3.16)  indeed  has  has  index  one  and  that  the  exact 
solution  is  x{t)  =  ((!  —  #)  e* -\-t^ ,e^  —  t^)^ .  LTVXE  performs  satisfactorily  for  all  tolerances. 
For  example,  a  run  with  LTVXE  from  #  =  0  to  #  =  8.0  with  tolerance  10“®  used  801  steps 
and  produced  at  #  =  8.0  the  approximate  solution  x  =  (—20354.512,  2916.9580)^  which 
has  a  relative  error  of  9.531  x  10“®  imder  the  maximum  norm. 
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